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We study discontinuous percolation transitions (PT) in the diffusion-limited cluster aggregation 
model of the sol-gel transition as an example of real physical systems, in which the number of 
aggregation events is regarded as the number of bonds occupied in the system. When particles are 
Brownian, in which cluster velocity depends on cluster size as ~ with r\ = —0.5, a larger cluster 
has less probability to collide with other clusters because of its smaller mobility. Thus, the cluster 
is effectively more suppressed in growth of its size. Then the giant cluster size increases drastically 
by merging those suppressed clusters near the percolation threshold, exhibiting a discontinuous PT. 
We also study the tricritical behavior by controlling the parameter 77, and the tricritical point is 
determined by introducing an asymmetric Smoluchowski equation. 

PACS numbers: 02.50.Ey,64.60.ah,89.75.Hc 



Percolation, a stochastic model for the formation of 
macroscopic-scale spanning clusters, has received consid- 
erable attention in statistical physics for a long time as a 
model for metal-insulator transitions, sol-gel transitions, 
epidemic spreading, fracture and so on [T]. When a con- 
trol variable, which is the occupation probability of a 
conducting bond between two vertices, is increased, a 
long-range spanning cluster emerges at a critical thresh- 
old "Pc- Such a percolation transition (PT) is convention- 
ally continuous. The discovery of a discontinuous PT has 
therefore been a long-standing issue in statistical physics. 
In this circumstances, a recently introduced stochastic 
model (Ij for the explosive PT has attracted considerable 
attention in a short time period. This stochastic model is 
a simple modification of the classical Erdos-Renyi (ER) 
random graph model [3], which contains suppression ef- 
fect in growth of cluster sizes. Such an explosive PT 
behavior has also been observed in other recently intro- 
duced stochastic toy models [iVfTU]. Even though it has 
been controversial if such an explosive PT is indeed dis- 
continuous in the thermodynamic limit pTMl3j . the ex- 
plosive percolation model has opened a new avenue for 
the study of discontinuous PT in non-equilibrium sys- 
tems. In this Letter, we are interested in what a physical 
system in real world a discontinuous PT can be observed. 

In this Letter, we consider the diffusion-limited cluster 
aggregation (DLCA) model [I4HT6] for the sol-gel tran- 
sition as a candidate in real-world systems of showing a 
discontinuous PT. This model was introduced a long time 
ago, and the dynamic cluster-size distribution were inten- 
sively studied for this model [T7lI21] . Here, this model is 
studied in the context of PT, which takes place as the 
number of cluster aggregation events increases. We also 
show that indeed the DLCA model exhibits a discontinu- 
ous PT. Furthermore, we generalized the DLCA model in 
which cluster velocity depends on cluster size as ~ s**. 
As the parameter 77 varies, there exists a tricritical point 
beyond which the PT becomes continuous with contin- 
uously varying exponents. We show that the general- 



ized DLCA model can be represented via an asymmetric 
Smoluchowski equation, by which the tricritical point can 
be determined. 

The DLCA model is simulated in the following way: 
Initially, single particles are placed randomly in L x L 
square lattices. Simulations start from N mono particles. 
The density of the particles is given as p = N/L'^. The 
system size L is controllable, while the density remains 
fixed in the simulations. Here we consider the case that 
particles are Brownian, so that velocity of a cluster is 
inversely proportional to the square root of its size (22) . 
To implement, we perform simulations as follows |16| : 
At each time step, (i) a s-sized cluster is selected with 
the probability q = s~°-^/(^^ iVsS~°-^), and is moved 
to the nearest neighbor. When two distinct clusters are 
placed at the nearest-neighbor positions, these clusters 
are regarded as being merged, forming a larger cluster, 
(ii) The time is advanced by 6t = l/(X]s -^sS~°'^), where 
Ns is the number of s-sized clusters. We iterate the steps 
(i) and (ii) until the giant cluster size is N. Later, we 
will consider a more general case in which the velocity is 
proportional to s'^ [TO] . 

When we study a PT problem of networks, control pa- 
rameter is the number of edges added to the system per 
the total number of nodes. Following this convention, we 
introduce a new variable p, which is defined as the num- 
ber of cluster aggregation events per the total particle 
number. Whenever two clusters merge, p is increased by 
6p — 1/N. Since N — 1 aggregation events occur during 
all aggregation processes, the aggregation event stops at 
Pf ~ 1 ~ 1/N. p depends on i in a nonlinear way as 
shown in the inset of Fig. [l|b) . 

In the original study of the DLCA model, the PT was 
not studied because the giant cluster size G{t) increases 
monotonically as the time t increases (Fig. [ij^a)). How- 
ever, we show here that when the giant cluster size G 
is traced as a function of p, it drastically increases, ex- 
hibiting a discontinuous PT, as shown in Fig. [l|b). This 
different behavior is caused by the nonlinear relationship 
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FIG. 1. (color online) (a) Plot of the giant cluster size G 
versus t when particles are Brownian. The giant cluster grows 
continuously from t = Q. (b) Plot of G versus p. G grows 
drastically near pf = 1 — 1/N. Inset: Plot of the relationship 
between t and p. t increases drastically but in a power-law 
manner as p approaches p / . Simulations are carried out with 
= 8,000 mono particles at t = on 400 x 400 square 
lattices. 



between t and p shown in the inset of Fig. [ijb). When 
p is small, t increases almost linearly with respect to p. 
However, as p approaches to Pf, t increases drastically in 
a power-law manner. As a result, G exhibits a discontin- 
uous PT. 

To check if the PT is indeed discontinuous, we use the 
finite-size scaling theory recently proposed for studying 
the explosive PT !5|. For the first, we measure Gjv(p) 
for different N under the condition that the particle den- 
sity p remains fixed. The results are shown in Fig. [2j^a) 
and (b). For a given iV, we pick up a Pc{N) at which 
the increasing rate dGN/dp is maximum. Then, the p- 
intcrccpt of the tangent of Gm{,p) at Pc, denoted by pd, 
is determined as 



Pd^Pc 



/ dGNjp) \ 
V dv T),/ 



dp 



Gn{Pc)- 



(1) 



Then, p^ also depends on N. We find that dGN{p)/dp 
at Pc increases in a power-law manner as ~ jyo.se 
(Fig.[2jc)), indicating that the giant cluster size increases 
more drastically as N increases. Thus, the transition 
is indeed discontinuous. Near pd{N), p is rescaled as 
p — {p — pd)dG m{Pc) I dp, which is then TV-independent. 
The giant cluster is then plotted as a function of p. In- 
deed, we can see that the curves of the giant cluster 
sizes for different N collapse well onto a single curve 
(Fig. [2jd)). Thus, the order parameter of the PT is writ- 
ten in the scaling form, 

G{p)^N-PI'h{{p-pd)N'"), (2) 
where fo{x) is a scaling function, and (3 — Q and 1/P w 
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FIG. 2. (color online) (a) Plot of Gat versus p for differ- 
ent particle numbers, TV = 50,000 (□), N = 98,000 (o), 
N = 162,000 (A), N = 242,000 (v), and TV = 338,0000 (o) 
at a fixed density p = 0.05. Gn begins to increase drastically 
near p « 1 and thus the data are not distinguishable for differ- 
ent sizes in the region far smaller than p — 1. (b) Zoom-in plot 
of Gjv versus p near pc- As the system size grows, the giant 
cluster grows more drastically, (c) Plot of dGN(p) / dp calcu- 
lated at Pc{N) versus TV. The slope dGN(pc) / dp increases as 
jY0.86±o.02^ indicating that it diverges in the thermodynamic 
limit, (d) Plot of G versus p = {p — pd)dG n(Pc) / dp for dif- 
ferent TV. Pd{N) is the p- intercept of the tangent of the curve 
Gn{p) at pe. (e) Plot of Xi/TV" ''^ versus p. (f) Plot of X2/N 
versus p. For (d), (e), and (f), the data collapse well onto a 
single curve. 



0.86 ± 0.02. We remark that this finite-size scaling form 
differs from the conventional one used in the continuous 
transition in the aspect that pd depends on the parti- 
cle number TV, which does in turn the system size L. 
Whereas, in the conventional scaling form used for a con- 
tinuous PT, Pd is replaced by P(.{oo), i.e., the critical 
point in the thermodynamic limit, which is independent 
of TV. For comparison, the modified ER models PH - 
I13| which were claimed to exhibit continous PTs show 
/3/z/ > even though their values are extremely small. 

We examine the behavior of the susceptibility. The 
susceptibility is defined in two ways. The first is the 
mean cluster size Xi{p) = Yl's '^^"■^(p)/ sns{ji), which 
exhibits a peak at pci- This quantity approaches pf de- 
fined earlier, as TV increases. The susceptibility at Pci 
increases with TV as Xi{Pci{N)) - 7vO-95±o.oi^ rpj^^^g^ 
Xi{p) is written in the scaling form Xiip) ^ N'''^/'^ fi{{p— 
Pd)N^^'^), where fi{x) is another scaling function and 
7i/p 0.95 ± 0.01. The scaling behavior is confirmed 
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numerically in Fig. [2]je). The other susceptibility is the 
fluctuation of the giant component sizes, i.e., X2{p) = 
N ^/ {G\{p)) — {Gn{p)Y- This quantity exhibits a peak 
at Pc2- We find that X2(Pc2{N)) - N. Thus, = 1 
and X2{p) N f2{{p — Pd)N^/^) with a scaling function 
/2. The scaling behaviors are also confirmed numerically 
in Fig.gf). 

The cluster aggregation process may be described via 
an asymmetric Smoluchowski equation: 



dp 



E 



Uskg 



C{p)C'{p) 



c{p) c'ipY 



(3) 



where = Ng/N is the concentration of s-sized clusters, 
which depends on p, and kik'j /{CC) is a collision kernel, 
where C{p) = Y^^kiTu and C'{p) = ^ik[ni. ki/C and 
k'j/C are the probabilities for i- and j-sized clusters to 
merge, in which the prime denotes mobile clusters and 
the other does immobile clusters in simulations. Their 
kernels are different below. 

We measure size-dependent behaviors of fc^/C and 
k'g/C numerically. Even though it is not manifest that 
ks/C and k'^/C follows a power law for the Brownian 
particle case (see Fig. |4]^a)), we roughly estimate that 
ki/C - and k'^/C - j-^ '^. The relation between 
these two exponents is mentioned later. With these ex- 
ponent values, we solve the asymmetric Smoluchowski 
equation numerically and obtain that /? = 0, = 1, 
71/1/ = 1 and = 1. These obtained numerical val- 
ues indicate more clearly that the transition is indeed 
discontinuous. 

We now consider a more general case in which the ve- 
locity is given as Vs oc [19], where s is cluster size. To 
implement this cluster is picked up with a proba- 

bility proportional to . The other rules in the numerical 
simulations remain the same. When clusters merge, the 
variable p is advanced by regardless of the cluster 

size s. The time is by 5t = ^/{J2b -^sS^)- Intuitively, 
when 77 is small or negative, fewer large-sized clusters are 
selected, so their growth is suppressed. Medium-sized 
clusters are abundant even close to pj, and they merge 
suddenly. In this discontinuous PT takes place. 

In contrast, when 77 is positively large, more large-sized 
clusters are selected and they have more chance of col- 
liding with other clusters, and merge to a bigger cluster. 
Thus, they can grow faster than smaller clusters can, 
so the giant cluster grows continuously and the PT is 
continuous. Thus we expect that there is a tricritical 
point across which the transition type is changed. We 
show snapshots of the system for different values of 77 at 
p = 0.99 in Fig. [3j which support the above argument. 

It has been roughly argued that the collision kernel in 
the Smoluchowski equation may be related to the perime- 
ter of a cluster as ki ^ i^-^l'^s [53]. However, when 
the selection probability is taken into consideration, the 
collision kernel of mobile clusters may be modified as 
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FIG. 3. (color online) Snapshots of the system for various 
values of ry at p = 0.99. The velocity of each cluster is given 
as i;^ oc s'', where s is cluster size. Numerical simulations 
are carried out for A?" = 8, 000 particles on L x L = 400 x 400 
square lattices. Since p is fixed, the number of clusters for each 
case is equal. The giant cluster is represented in a different 
color (gray, orange). The cluster-size distribution becomes 
more heterogeneous as rj increases. 
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FIG. 4. (color online) (a) Numerical estimations of k[/C' at 
Pd for -q = -0.5 (□), r; = (o), 77 = 0.4 (A) and r? = 0.8 
(V). Simulations are carried out for = 8,000 and L = 
400. Slopes of the guidelines are -0.23 ±0.02 for (r? = -0.5), 
0.32 ± 0.04 (»7 = 0), 0.62 ± 0.01 (7? = 0.4), and 0.88 ± 0.01 
{rj = 0.8). (b) Numerical estimations of the ki/C at pd- The 
same symbols are used as (a). System size is TV = 8000 and 
L = 400. Slopes of the guidehnes are 0.35 ± 0.04 (77 = -0.5), 
0.32 ± 0.04 (?7 = 0), 0.2 ± 0.01 (ry = 0.4), and 0.05 ± 0.005 
(r, = 0.8). 



k[ ~ i**^^ 1/''/. For the Brownian case with 7/ = —0.5, 
by using the fractal dimension df = 1.4 — 1.5, the mea- 

' ' ,-0.3 



sured values H 



for mobile clusters and fc, 



for immobile clusters (Fig. |4|^a) and (b)) are reasonable. 
When rj is sufficiently large such as 77 > 0.8, the largest 
cluster grow by merging small-sized clusters. In such a 
merging process, small-sized clusters can penetrate into 
the interior of the giant cluster, and then the merging 
probability can be independent of cluster size (Fig.|4j^b)). 
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FIG. 5. (color online) Plot of Gjv versus p for the Smolu- 
chowski equation with fci ~ 1 and k'j ~ j" ** (a), and fci ~ 1 
and k'j ~ j^'^ (b). As the system size grows, Gn grows more 
drastically and Pc{N) approaches one in (a), but it decreases 
to Pc(oo) > in (b). For (b), the transition is continuous. 
Similar plot for the DLCA with rj = 0.8 (c) and rj = 1.5 
(d). Similar behaviors are shown. Simulations are performed 
for N/W^ = 8,32,72,128 and 200. Inset of (d): Finite size 
scaling of Gn{p) for different system sizes with l/i/ = 0.28, 
P/v = 0.5 and pc(oo) — 0.42. These numerical values depend 
on rj. Thus, for rj > rjc, the critical points for different rj form 
a critical line. 



Thus, k[ ^ i^' and ki ^ constant. 

We integrate the Smoluchowski equation with using 
k'JC ~ and h/C ~ 1. We find that the transition 
behavior changes across 7]^ « 1.3 — 1.4. When r] < rjc 
iv > Vcji as the system size increases, the percolation 
threshold increases to one (decreases to a finite perco- 
lation threshold). Those behaviors can be observed in 
Fig.jsjja) and (b). Moreover, for the latter case, the tran- 
sition turns out to be continuous. Therefore we con- 
clude that there exists a tricritical point which locates 
at rjc- Similar behavior is observed for the DLCA model 
(Fig. [sjjc) and (d)). We also check the cluster size distri- 
bution at the percolation threshold. Indeed, the distri- 
bution obeys a power law Us ^ s^"^. Also, the exponent 
satisfies t < (>)2 for 77 < {>)ric ~ 1-3 - 1.4. The ex- 
ponent r = 2 is marginal between the continuous and 
the discontinuous PT, which was proven analytically in 
the cluster aggregation model Therefore our numer- 
ical result confirm that the tricritical point locates near 
77e«1.3-1.4. 

In summary, we have studied the DLCA model as an 
example of real-world systems exhibiting discontinuous 
PT. The velocity of an s-sized cluster is given in a gen- 
eral form Vg ^ s^. When rj < {>)r]c, the PT is discontin- 
uous (continuous), where a tricritical point rjc is roughly 
estimated to be 1.3 — 1.4. Since the case 77 = —0.5 cor- 
responds to the Brownian particle motion in real-world 
systems, we can say that a discontinuous PT can take 



place in real-world non-equilibrium systems. We also in- 
troduced and studied an asymmetric Smoluchowski equa- 
tion, and determined the tricritical point from the fact 
that the cluster size distribution follows a power-law be- 
havior with exponent —2 at the tricritical point. We 
finally remark that the discovery of the explosive PT in 
the DLCA model was made by tracing the giant clus- 
ter size as a function of p instead of indicating that 
a discontinuous PT may be explored as a function of a 
unconventional parameter. 
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